RESIDUAL DISTRIBUTION SCHEMES FOR CONSERVATION LAWS VIA ADAPTIVE 

QUADRATURE 

TIMOTHY BARTH* AND REMI ABGRALL* 

Abstract. This paper considers a family of nonconservative numerical discretizations for conservation laws which retains 
the correct weak solution behavior in the limit of mesh refinement whenever sufficient order numerical quadrature is used. Our 
analysis of 2-D discretizations in nonconservative form follows the 1-D analysis of Hou and Le Floch [14]. For a specific family of 
nonconservative discretizations, it is shown under mild assumptions that the error arising from nonconservation is strictly smaller 
than the discretization error in the scheme. In the limit of mesh refinement under the same assumptions, solutions are shown 
to satisfy an entropy inequality. Using results from this analysis, a variant of the W N” (Narrow) residual distribution scheme of 
van der Weide and Deconinck [31] is developed for first-order systems of conservation laws. The modified form of the N-scheme 
supplants the usual exact single-state mean-value linearization of flux divergence, typically used for the Euler equations of 
gasdynamics, by an equivalent integral form on simplex interiors. This integral form is then numerically approximated using 
an adaptive quadrature procedure. This renders the scheme nonconservative in the sense described earlier so that correct 
weak solutions are still obtained in the limit of mesh refinement. Consequently, we then show that the modified form of the 
N-scheme can be easily applied to general (non-simplicial) element shapes and general systems of first-order conservation laws 
equipped with an entropy inequality where exact mean-value linearization of the flux divergence is not readily obtained, e.g. 
magnetohydrodynamics, the Euler equations with certain forms of chemistry, etc. Numerical examples of subsonic, transonic 
and supersonic flows containing discontinuities together with multi-level mesh refinement are provided to verify the analysis. 
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1. Motivations. Discrete conservation has become a standard design criteria in the development of 
numerical discretization techniques for conservation laws that admit discontinuous solutions. From the Lax- 
Wendroff theorem [21], the ingredients of consistency, stability, and discrete conservation yield convergent 
approximations of conservation laws in divergence form for both smooth and discontinuous solutions that 
are valid weak solutions in the sense of distribution theory. Even so, the development of stabilized numerical 
discretizations often also utilizes the quasilinear form (a.k.a. nonconservative form) of the conservation law 
system to approximate simple or plane wave solutions for use in upwind stabilization mechanisms. As we 
will illustrate later, the use of quasilinear forms is often at odds with the requirement of discrete conservation 
unless specialized mean-value linearized variants of the discrete quasilinear form are used. As a practical 
matter, obtaining simple expressions for these mean-value linearizations in closed form is often extremely 
complicated or even impossible. In addressing this difficulty, our goal is to develop a general framework that 
avoids these complications while still insuring that valid weak solutions of the conservation law system are 
obtained in the limit of mesh refinement. 

As a motivating example, consider the scalar Cauchy problem in one space dimension and time 

f u tt -F (/(u)),* = 0 for (x, t) €Rx R + n.l) 

\ u(x,0) = tx 0 (x) 

with uGl and f(u) : R R. In this equation u 0 (x) is assumed to be periodic or compactly supported 
data. Let Ax j+ i/ 2 = x J+ 1 — xj denote a general nonuniform partitioning of space so that Uj represents the 
numerical approximation u(xj,t). Next, consider the prototype conservative semi-discrete scheme 

dUj ftj+l /2 " hj- 1/2 _ Q ( L2 ) 

dt Axj 

with hj±i / 2 the numerical flux. This prototype scheme is conservative in space due to the mutual telescoping 
of numerical fluxes. A first order accurate upwind scheme is easily obtained via the flux function 

h j+i/2{uj,u j+1 ) = i ( }{uj ) + f(u j+ i)) - ^|a| j+ i/2 (tij+1 - Uj) (1.3) 
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with dj+i / 2 an approximation of the flux Jacobian df/du at Xj + 1 / 2 - Observe that whenever the exact 
mean-value linearizations are used, e.g. 


f(uj+ 1) - I(Uj) = (o)j + 1/2 (ttj+l - u j) 

so that Uj±i /2 = {a) J±1 / 2 , the first order upwind scheme can be written equivalently as 




~ u i j. /„\+ ~ Uj ~ L = 0 . 


(1.4) 


(1.5) 


Note that this discretization is nonconservative in space unless the exact mean- value linearization (1.4) 
is used. Nonconservative schemes of this form axe known to converge to incorrect weak solutions. More 
precisely, Hou and Le Floch [14] have shown (in 1-D) that if the nonconservative scheme (1.5) converges, it 
converges to a solution of 

+ (/(u)), x = A* 

where /i is Borel measure source term that is expected to be zero in the regions where u is smooth and 
concentrated where u is not smooth. The construction of an exact mean- value linearization is readily 
accomplished in 1-D by the general path integration 


/*tl B [UB 

/(ufl)-/(tM)= / df= a(u)du 
Ju A JUA 


-l 


U A 

■Cb 


U A 

a(«(0) ^ d Z - 


( 1 . 6 ) 


Without loss of generality, one can restrict u(f) to the space of polynomials, P*. A particularly convenient 
choice are Pi linear polynomials since 


du 

f{u B ) ~ f{u A ) = J f a(u( 0 ) 

= f a(u(0) d Z 

JiA 




«(0€Pi 


( UB - UA \ 

\eB-u) 


(1.7) 


so that the following mean-value speed is obtained 


(a){u A ,u B ) = 1 — [ a(u(f)) d£ 

” <i A J^A 


( 1 . 8 ) 


tiU)€Pi 


In Harten, Lax, and van Leer [13] this expression is interpreted as an integration in state space parameterized 
along the line 7ru(0 = im + f i u B — tM), £ £ [0, 1]. 


(a)(u A ,u B ) = [ a(nu{ 0 ) d Z * 
Jo 


When the locations A and B are not coincident, one can equivalently interpret this as an integration in 
physical space assuming the Pi Lagrange interpolant 


iru(x) ~ u A H — (uj 3 — u A ), x G [x Ay x B ] 


x B ~ x A 


so that £ = x and 


x B 


1 f XB 

/ a(7nx(x)) dx . 

~ J X A 


(a)(u A ,u B ) = 


(1.9) 
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This latter interpretation is useful since it generalizes the mean-value contruction to simplices and more 
arbitrary regions. Next consider an approximation of Eqn. (1.9) using NQ- point numerical quadrature 

NQ 

{q)(ua,v>b) = Yj ui a(nu(qi)) + Rnq + i ( 1 - 10 ) 

i=i 

where are quadrature weights, qi quadrature positions, and Rnq+ i is the numerical remainder term. This 
renders the scheme (1.5) nonconservative in space. In later sections, we derive (under suitable assumptions) 
the same result as Hou and LeFloch and are able to characterize more precisely the Borel measure fi. In 
particular, the dependency of fj. with respect to the number of quadrature points is given. If an adequate 
number of quadrature points is taken, the error terms due to nonconservation are shown to be comparable 
or smaller than the discretization error of the scheme. In addition, a discrete entropy inequality is formally 
obtained in the limit of mesh refinement. From the practical point of view, these results are important with 
the following consequences 

• Exact mean- value linearization is no longer needed. This is useful when solving systems of conser- 
vation laws for which exact mean-value linearizations are not known in closed form, e.g. magneto- 
hydrodynamics, Euler equations with certain forms of chemistry, etc. 

• General finite element shapes are permitted, e.g. tetrahedra, hexahedra, prisms, pyramids. Previous 
exact mean-value linearizations in closed form have been restricted exclusively to simplex shapes. 

The new nonconservative formulation suggests an adaptative strategy, whereby the number of quadrature 
points depends on the local smoothness of the numerical solution. This strategy is undertaken in Sect. 3. 

2. Background. In this section, we briefly review a number of well known constructions and analytical 
results that we utilize later in the development and analysis of our nonconservative formulations. 

2.1. Conservation Laws and Symmetric Hyperbolic Forms. Consider the Cauchy problem for 
m coupled conservation laws in d space dimensions and time 

I w,e+ E r(w) iX . = 0 for (x,t) £ R d x R+ (2 .i) 

[ w(z, 0) - W 0 (x) 

where w £ denotes the vector of conserved variables and f(w) : IR m IR mX£ * a flux vector. In addition, 
Eqn. (2.1) is assumed to be equipped with a convex entropy extension so that the additional scalar inequality 
holds 

<° (22) 

i— 1 

with H( w) : ]R m M the convex entropy function and G(w) : W nxd H > the entropy flux vector for the 

system. Solutions of Eqn. (2.1) satisfying (2.2) are generally of two types [22]: 

• (Classical Solutions) Smooth solutions satisfying the quasilinear form of Eqn. (2.1) 

d 

w it + Y A,(w) w |X . = 0, Ai{ w) = P w . (2.3) 

i=l 

As part of the symmetrization theory for first-order conservation laws developed by Godunov [11], 
Mock [23] and others, it is known that the existence of a convex entropy extension insures that the 
quasilinear form (2.3) is symmetrizable via a change of variables w v where v = H^, £ IR m 
denotes the so-called entropy variables for the system. As consequences of symmetrization theory, 
performing the change of variables 

d 

Ao v < + AiV, Xi = 0 , 


(2.4) 


yields the matrix A 0 = w, v = (H, w , w ) -1 symmetric positive definite and the matrices A, _ f ‘ v - 
AiAo symmetric. For brevity, the functional dependency of the matrices -4, and A t has been omit- 
ted. Motivated by the energy analysis given in subsequent sections, we assume the basic solu- 
tion unknowns are the entropy v-variables so that the shorthand notation .4,(w) should be inter- 
preted as Aj(w(v)). It is useful for later developments to define the real- valued matrix combination 
A( w,w) = WjAi(w), oj G R d and similarly the symmetric matrix i(w,w) = WiA(w). Observe 
that symmetry of A( w,w) implies that A(w, w) has m real eigenvalues, A! < A 2 < • • ■ < A m and m 
real eigenvectors rfc(w,ui) 6 R m satisfying the standard eigenvalue problem 

A(-w,u>) rjt(w, w) = Ajfc(w,w) r*(w,w), k = 1,2 ,...,m 


since the identity 

i" 1/2 Al(w,w)iJ /2 = i" 1/2 i(w,w)io 1/2 
shows that ,4(w, ui) is similar to a real- valued symmetric matrix. 

Our keen interest in the quasilinear form (2.3) comes from its use in the construction of upwind 
discretizations such as variants of Godunov’s method [10] utilizing approximate Riemann solvers 
[32, 26] and the multi-dimensional fluctuation splitting scheme described in the following section. 
Specifically, the quasilinear form (2.3) admits nonlinear simple wave solutions of the following form 
for a given direction vector cj: 

m 

w(x, t) — * x — Ajt(W fc ,cj) t)) (2-5) 

*=i 

where W k {o,u) 6 satisfies the differential relation 

~ =T k (W k (cr,u),u) ( 2 - 6 ) 

do 

for the self-similar real-valued parameter a. In Eqn. (2.5), a* G R are expansion coefficients to 
be determined by matching initial data. When the matrix .4(w) is assumed locally independent of 
w, then local plane wave solutions are obtained. Historically, mean- value linearized variants of the 
quasilinear form (2.3) have been used in 1-D to construct approximate Riemann solutions [26] for 
eventual use in upwind discretizations. In Sect. 2.2, we consider a multi-dimensional upwinding 
strategy which also uses plane wave information originating from a mean-value linearized form of 
Eqn. (2.3). 

• (Discontinuous Solutions) Weak solutions of the divergence form (2.1) satisfying a jump condition 
on space-time hypersurfaces, S, with space-time normal vector n = («t,n T ) 

n t [w]t + ^2 n *[^]- = 0 

i=l 

with [arg((x,t))]± = lim ti0 (arg((x,t) 5 +en) - arg((x, t)s - en)). In Sect. 3, a Lax-Wendroff-hke 
theorem is presented which addresses the convergence to weak solutions of a family nonconservative 
discretizations using approximate mean-value linearization. 

Note that in the remainder of the paper, the notation || • || will denote a pointwise norm over m variables 
unless otherwise indicated. When the argument is dime nsionally comparable with the v-variables, the 
natural norm is not the standard Euclidian norm ||x|| = VE'.li *1 but rather the dimensionally consistent 
matrix norm [15] 

||x|& 0 = x T io* 

where Aq is the inverse of the Hessian matrix of the entropy, -lo = (ff,w,w) 


( 2 . 8 ) 
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2.2. The Residual Distribution Scheme on Simplicies. In the remaining sections, we assume a 
triangulation T h in R d of a polygonal spatial domain f! composed of nonoverlapping simplices T iy Q = UT,, 
Ti D Tj = 0, i ^ j. A simplex T in R d is uniquely described by d + 1 vertices T{M U M 2 , • • • , M d+l ). For 
purposes of analysis, the triangulation is assumed to be shape regular with maximum simplex diameter h. 
From the triangulation Th, the geometric dual tessellation Ch is constructed by connecting gravity centers of 
the simplices and the mid-points of the edges as shown in Fig 2.1. In this figure, the dual cell Ci surrounds the 
triangulation vertex Mi. We also define piecewise linear and piecewise constant spaces on the tessellations 



Fig. 2.1. Dual cell C{ associated with triangulation vertex M x in R 2 


Th and Ch respectively 

v„ = {v fc ;v h € C°(R d ) m ,v h|T € (?>i) m ,V T 6 %} 


x h = {vhjvhic e CP 0 ) m ,v c e c h ] . 

Let V, € R Tn denote the nodal degrees of freedom located at Mi which uniquely describes v/, in both 
spaces V h and X h . For example, if N,(x) denotes the standard piecewise linear basis function for triangulation 
Th such that Ni(Mj) = <5 tJ , then for v h € Vh 

v*(x)= £ Ni(x)Vi . 

MieTh 

Similarly, if Xi( x ) denotes the characteristic function for the dual cell Ci € Ch, 



x £ Ci 

x&Ci 


then for £ Xh 

v 'h( x ) = 5Z • 

Mi&n 

Finally, for brevity of notation, we shall write w h = w(vh) and W , = w(V,) to denote the corresponding 
conserved variable forms. 

Using these definitions, we can state the simplest prototype residual distribution scheme (explicit in 
time) used in discretizing (2.1). 

Residual Distribution Scheme: For all Mi £ Th and n > 0 
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where € lR m represents a discretization of the negated time evolution term integrated in the simplex T 




-b 


W h),t dx 


(2.10) 


with an as yet unspecified sum decomposition of among the d + 1 vertices of the simplex T in R. d 

= $l,T + ' ' ' + $d+l,T • (2-H) 


In the special case of Eqn. (2.1), Is expressed equivalently in terms of the spatial flux divergence 

+ 1" $d+l,T := J ^53 dx . (2.12) 

The residual distribution scheme encompasses a number of well known weighted residual methods that have 
residual decompositions which reduce to following form for P\ linear elements: 

* i,T = (dTT + ft.13) 

where r € R mxm is an nonsingular matrix such that Yli=i T = °* Some examples of weighted residual 
methods for solving (2.1) include 

• The streamline diffusion method of Johnson and coworkers [17, 18]. 

• The streamline upwind Petrov-Galerkin (SUPG) and Galerkin least-squares finite element methods 
of Hughes and coworkers [15, 16]. 

• The cell vertex finite volume methods of Ni [25] and Morton et al. [24, 7] 

The residual distribution formula also describes a family of monotone and positive coefficient schemes for 
scalar conservation laws due to Roe [27, 28] and Deconinck [8] and the system extension due to van der 
Weide and Deconinck [31] described later in Sect. 5. Fundamental to these residual distribution schemes is 
the mean-value linearization of the flux divergence formula (2.12). 


^^(w),,^ dx = ^2(A)i ( 
,i= i / «=i Jr 


w 


dx 


(2.14) 


To facilitate this calculation, we follow the 1-D example of Sect. 1 by introducing an auxilliary mapping 
z(v) : lR m *-¥ lR m and restricting z to the space of piecewise linear Lagrange interpolants denoted by tt^z. 


J (y: f 1 (w(7T/ 1 z) |X .^ dx = y (i4)j J w(7T h z), Ij dx (215) 

For the Euler equations of gas dynamics, the choice z = ( y/p , y/pV , y/pHt) T with p the fluid density, V the 
fluid velocity, and Ht the fluid total enthalpy yields closed form expressions for the mean- value linearization 
of the flux divergence [30]. A striking property of this linearization is that the linearized system 


w,* + YMh".** = 0 ( 2 * 16 ) 

t=i 

is hyperbolic. Unfortunately, no other z variable is known to give simple and closed form formulas leading to 
a hyperbolic linearized system. In addition, this approach is limited to simplices, while in many applications 
hexahedral brick meshes would be desirable for accuracy reasons. 

From a theoretical and practical point-of-view, there is motivation to work directly with the entropy 
variables since the corresponding mean-value linearized form 




v, x < dx 


(2.17) 
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would necessarily produce a hyperbolic linearized system due to the symmetry of (A)j. Using symmetric 
forms, we also show in subsequent analysis the satisfaction of an entropy inequality in the limit of mesh 
refinement. Our general strategy is to utilize a piecewise linear representation of the entropy variables 
themselves so that v/* € V* and 


L(t« w(v*))^) dx = \ dx = l T l ( 218 ) 

Jt \ i=l / «=1 JT »=1 


with 

{A)i = E f r Mvh) dx . (2.19) 

Following the 1-D motivational example of Sect. 1, Eqn. (2.19) is approximated by quadrature formula so 
that componentwise 


NQ 

(A),- = + Rnq+ i • (2.20) 

i=i 

Since \h\r € Pi (T), Eqn. (2.18) has used the fact that the gradient component are constant within a 
simplex. Consequently, the quadrature formula used in (2.20) 

NQ 

/ H(x)dx = \T\J2u>iH( qi )+0{h k + l ) (2.21) 

'* T <=i 

should at least be exact for H(x) € Pfc(T) and k > 1. In addition, the 0(/i* +1 ) error is assumed to have the 
following behavior for use in later analysis: there exists C independent of the simplex T such that 

0(h k+l ) < C(%) j T \\D k+1 H(x)\\dx (2.22) 

where h is the maximum diameter of the T, C(Th) is a geometrical parameter that only depends on 7^, and 

DkR{x) = {^' |a| = *} • 

Note that the use of numerical quadrature permits generalization of the techniques to non-simplicial elements, 
e.g. brick elements (Q), using the form 

NQ / d \ 

*q = IQI ( 'Y^Ai[y h {.qi))(yh{qi)),x i I +Rnq+i • ( 2 - 23 ) 

1=1 \i=l / 

Hence, we are interested in residual distributive schemes that fullfil the approximate conservation relation 

NQ / d \ 

*1,T + ' • ' + ^+1,T = l^| E ^( V »(«))( V ‘(«)).'. • ( 2 ‘ 24 ) 

|=1 \t=l / 

In Sect. 5, a particular residual scheme known as the N-scheme is considered [30] as generalized to 
systems of conservation laws by van der Weide and Deconinck [31]. This system N-scheme assumes an 
exact mean- value linearization via the parameter vector. We then propose a variant of the system N-scheme 
which utilizes a piecewise linear space consisting of the entropy variables and approximates the mean-value 
linearization via quadrature. Analyzing this new scheme for systems of conservation laws, we show that in 
the limit of mesh refinement that numerical solutions satisfy an entropy inequality. We then show a similar 
result for the system N-scheme when the linearization is approximated via quadrature. 
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3. Weak-* Convergence, a Lax-Wendroff Result. Consider the numerical scheme (2.9). The 
nodal variables WJ 1 are assumed to map uniquely via v(w) and w(v) to and from Vj* which are the degrees 
of freedom in the spaces V/, and X h at time t n = nAi, n € [0,iV]. In addition, the as yet unspecified residual 
decomposition and V" are assumed to satisfy the following conditions : 

ASSUMPTION 1 (HI). Let Tn be a shape regular triangulation. For C G R and any fixed n , there exists 
C f (C) G R which depends on the triangulation Th such that V v£ G Vh and ||vJ|| i0 o(Rd)m < C 

ll*”r|| <C'h d ~ l ||Vy-V?||, VTeT h andVMitT . (3.1) 

Mj GT 


This is a continuity assumption on the residual decomposition in a simplex T in terms of the local nodal 
values ofVfyMi G T. In particular , whenever v£ is constant in T we then require that $£ 7 * = 0. 

Assumption 2 (H2). For all v£ e Vh and fixed n 

d+l NQ / d \ 

*?• = E = PI E "M E MW*)) (vj(*))^. ) - « 6 T (3.2) 

i=l 1 = 1 \t=l / 


where A{ = f* v and NQ denotes the number of quadrature points. In addition, the quadrature error in the 
flux divergence calculation is assumed to be of the following form for all n and a given integer k > 1 


*r-j T ±r 

Jl 1=1 


? x ,(vh)dx 


< C(Th) 


h k+l 

(k+iy. 


D k+1 




(3.3) 


by using a sufficient number of quadrature points. 


Remark 1. 

(i) Observe that Vh G Vh is C° continuous, consequently for neighboring simplices sharing a common spatial 

edge, T jk = {x | dTj n dT k ^ 0}, 

pWWtlr, -nf = E fi KW)lr fc , * € T jk (3.4) 

1 = 1 1=1 

where is a directed normal on r^. 

(ii) For any constant C and fixed n, there exists C ( {C) such thatVv n G V^; ||vj^ | < C- Consequently , 

for shape regular T G Th 


*rii<x E ii v ? - v rn • 


(3.5) 


Af it Mj€T 


(iii) Lastly, for any sequence (v£)/> such that (v£) is bounded in L°°( R d x ) m independantly of h and N 
and converges in L / 2 0C (R <i x K + ) m to v, we have 


Hm o ||r( v h )-r( v )|| tL( R JxH+ ) m =°, i = 1,2, ...,d 


(3.6) 


Our first principle result is a generalization of the Lax-Wendroff theorem to residual distribution schemes 
for systems of conservation laws using numerical quadrature. Note that since the mapping H( w) is smooth 
and w h is bounded, assumptions on are equivalent to the same assumptions on defined by the nodal 
values V*. 

THEOREM 3.1. Consider an initial condition v 0 G L°°(R d ) m for time r > 0. Let Wj be the nodal 
approximation for all Mi G Th given by (2.9) from which Vj are obtain via V* = v(Wj). Assume that the 
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scheme satisfies assumptions (HI) and (H2) and that there exists a constant C that depends only on vq and 
functions v E L 2 ( R d x 1+ ) m and v* such that for v* E Vh 

supsup||v A (x,t)|| < C, lim ||v A - v|| L? (RrfxR+)m = 0 . (3.7) 

h x,t ,oeV 


Let Q = U T be a bounded domain of R d and t > 0 a bounded time. Assume that there exists a locally 
bounded, positive measure p such that ||Z?v/,|| tends to [i in the sense of distributions as h — ► 0. Then v(x, t) 

satisfies 



ip(x, 0)w(v o (x))dx 


< 


C(T h , f) 
(fc + 1)! 


(M,m> 


(3.8) 


where k is a integer as described in Assumption (H2) and C(Th, f) is a constant that depends on Th and 

ll^ +1 f,v||. 

This results applies when the limit is piecewise smooth, as it is in practical applications. The proof is 
inspired by [20] then [2]. 

4. Proof of Theorem 3.1. The proof of Theorem 3.1 appeals to a sequence of lemmas that are 
somewhat classical but axe tailored here specifically to residual distribution schemes and the use of numerical 
quadrature for element interior integrations. For simplicity, we assume an evolution to time r, an N integer 
multiple of At, i.e. r = N At, although the generalization to arbitrary bounded values of r is straightforward. 

Lemma 4.1. Let Q = U T denote a bounded domain of R d and t > 0 a bounded time . Further , let (vh)h 
denote a sequence such that v/j( . ,nAt) E for any n < N. Assume there exists a constant C independant 
of h and N and a function v E L 2 (Q x [0,r]) such that 

supsup||v h (x,l)|| < C, lim ||v h -v|| L? (Qx[0 ])m =0 . (4.1) 

h x,t n-»0 loc 

Under these assumptions , the following limits and bound are obtained 

1. lim*_ 0 El o At £vr eC \T\ Z Mi , Mj€T ||V? - V*|| = 0. 

2. lim/j-KO En=o ^ Hvrec 1^1 £m,,a/j€T ||^? — V"|| =0. 

3. lim A _ 0 h ||/?«Vfc|| L a ( Q X(0fT])m =0. 

4 . There exists C f independant of h and n such that /i||Z? a .v/ l ||£ap(Q X [o, r ])"» < C * . 

Proof To prove this Lemma, one needs only consider real-valued functions. For any simplex T and 
open time interval I n =]t n _i,t n [, two piecewise constant functions can be constructed which are useful in 
analysis, namely 


v/,(i,t)|Tx/’> = E XCinr(x)V" 

M iGT 


(4.2) 


and the shifted variant 


v k (*,0lTx/-= E ^nr(x)V" (0 (4.3) 

Mi€T 

where a(i) denotes a cyclic permutation of the index i and XCiHT is the characteristic function of C, O T 
with Ci the dual cell at node Ai*. This defines two functions on Q x [0, r] that are bounded independantly 
of h and N. Moreover, the following useful identity holds for these functions in a simplex T for arbitrary 

p > 0 


\T\ E liv?-vyir = (d+l) f \\v h -^h\\ p dx 


(4.4) 
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where the “d + 1” factor comes from the definition of dual cells in III d . Integrating in time 

Y At Y m y n v "- v "ii = r f iK-vhiidxdt . 

A W'T/- n i / TT^ Jo J(UT)eQ 


(4.5) 


n= 0 VTGQ 




The sequence (vh)h is bounded, therefore a function v 7 E I/°°(Q x [0,r]) m exists such that v* v 7 for the 
weak-* topology. From the previous assumptions, v/» -> v in Lj* oc which implies v = v 7 since Q x [0, r] is 
bounded and Cq^Q x [0, r]) is dense in L l ( Q x [0, r]). Similarly, there exists a function v E L°°(Q x [0, r]) m , 
such that v/j — > v in the weak-* topology. Our next task is to show that v = v and thus finally v = v = v 7 . 
To do so, let p(x , t) E Cq° ( R d xl + ), integrate p v/» in Q x [0, r] , and use the definition of the shifted function 
v 


[ [ ¥> v f> dxdt = Y At Y V <* [ VXCinTdxdt 

^ 0 n=0 VTGQA/.GT , ' T 

= f [ dxdt + V At V V t n f [ pxdnT dx- f ¥>Xc r -, m nT dx") (4. 

Vo Jq “ 0 ^ Ur Jr / 


6 ) 


where <7 x (i) denotes the inverse cyclic index permutation such that a(p l (i)) = i. Due to use of gravity 
centers and edge mid-points in the definition of the dual cells Ci, 


J XCiHT dx = J Xc v - i (i) dx = \Ci fl T|, i = 1, 2, . . . ,d+ 1 . 

Using the integral mean- value theorem, points X* and x 7 E T can be found such that 
Jr <P XCiHT dx = \Ci n T\ <p(xi), p xc„- 1 (0 dx = \CiC[ T\ i^(x') 
Since ||D<£>|| is bounded on Q x [0, T] and V* is bounded, 


[ [ pw h dxdt - f f 

Jo Jq Jo Jq 


ifVh dxdt 


< Ch 


(4.7) 


(4.8) 


(4.9) 


where C is independant of h and N . Hence in the limit v = v and finally v = v = v 7 . 

Let Vhy vtn and v' h denote scalar components of the respective vector- valued functions v/j, V/,, and v^. 
By the same technique, we see that the components (u£) and (vfj have the same weak-* limit. It will now 
be shown that this limit is v 2 . Once again appealing to the density of Cq°(Q x [0, r]) in L l (Q x [0 ,r]) and 
the fact that v\ is bounded independantly of h and N , we will take test functions p in C§°(Q x [0,r]). The 
function ip is bounded in Q x [0,r] and ► v in L 2 oc (Q x [0, r]), thus 


L 


CxJO.r] 


ip\v — t;h| 2 dxdt — ► 0, 


(4.10) 


(4.11) 


and consequently 

/ v 2 dxdt - 2 (pvvh dxdt + I p v\ dx dt -> 0 . 

Jqx[0,t] J Qx[0,t] J Qx[0,t] 

By the Cauchy- Schwarz inequality and <pv E L l (Q x [0,r]), the second term converges to 

f pv 2 dxdt . (4-12) 

J Qx[0,r] 

Hence, v % — t/ 2 in L°° weak-*. We are free to choose p = 1 combined with the limit -> u 2 in L°° weak-*, 
yielding 


L 


Q*(0,r] 


|57/i — u| 2 dxdt — > 0 


(4.13) 
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and finally 


/ \vh - Vh\ 2 dxdt -> 0 . (4. 14) 

J Qx[0,r] 

Interpreting this equation of the form (4.4) gives the asserted limit (1) of Lemma 4.1. The limit (2) of Lemma 
4.1 is then clear: Q x [0, r] is bounded thus L 2 (Q x [0, r]) C L l {Q x [0, r]). To prove limit (3) of Lemma 4.1, 
consider Fig. 4.1 which shows a simplex with inward pointing normals scaled by the edge length. In , the 
normal n* is the inward pointing vector pendicular to the ( d — l)-dimensional simplex facet opposite vertex 
Mi, i = 1, 2, . . . , d + 1 scaled by the measure of this facet so that Yli=i = 0* Using this notation 


3 



Fig. 4.1. Depiction of inward pointing normals, n* , for a simplex in R 2 . 


D x Vh\T — 


1 


d + 1 


d+1 


(d + l)\T\ 




3 = 1 


(d+\)\T\*z 




(4.15) 


Integrating in time and space 
N , , N 


E/ / ipxv^ii 2 dxdt^^At e m||(^ v h)ir|| 2 <c^E At E E i T ill v ”- v " 

n=0 JlnJ Q n = 0 VT€<2 n=0 VTeQ Mi,Mj€T 


because the the gradient is constant within a simplex and the triangulation is regular. Limit (3) of Lemma 
4.1 is then obtained from the application of limit (2). To obtain the bound (4) of Lemma 4.2, we again 
consider Eqn. 4.15 assuming bounded Vj 


h ||L > x v/ l ||/ <O o( ! 2 X [ 0ir ])m < 


h 


max 


(d + 1) |T| i<j<d+i 


| fp 


max V,* < 

l<j<d+l 


max 


(d + 1) |T| 


(4.16) 


which is bounded from above by a constant independent of h and N for shape regular triangulations. This 
concludes the proof of Lemma 4.1. □ 

Lemma 4.2. Let ip(x ,£) E Co(R d x R + ). With the assumptions of Theorem 3.1 , we have 


N 


E A * E l^^( M i.<n)(W" +l 

n— 0 Af , e Tk 


W") + [ ^£(x,t) -w h (x,t) dxdt + f <p{x,0) w 0 (a:) dx -> 0 

J R rf xR+ V* 

(4.17) 


when h -¥ 0. The proof is classical, see for example Kroner [19], p. 377. 

Lemma 4.3. If v/ l (x,t) € V/» satisfies the assumptions of Theorem 3.1 , then for any bounded Q and 
smooth ip(x, t), 


lim sup h k + l 
h 


E E [ Dz +1 (EA( v ft)( v *).-«) 

n— 0 VT6Q *' T \*=1 / 


dx 


<C(7i,f)(M,M> 


(4.18) 
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where 7r/i^(iT,t n ) is the midpoint value of the linearly interpolated <p function in simplex T for constant t n 
and C(Th,f) is a bound on ||D* +l f V (v/j)|| for bounded v^. 

Proof : Using the bound ||D* +1 f |V || < C{ 7i,f) together with the observation that D x v h is constant in a 
simplex T 


h k+l \\D k x +1 (^Ai(v h )(v h ), x J\\ = h k+l \\D k v +l (Y^Mv h ))((v h )' Xi ) k+2 \\ 

<ft* +1 ||^ +1 i(v fc )|| ||PxV h || fc+2 

<C(n,f, v )\\D x v h \\ . (4.19) 


Using again the fact that D x \^ is constant in a simplex, it follows that 

N f N J. 

hk+l Y, E *fcV»(*T,*n) / ||D* +1 .A(v k )|| ||Z?xV h ||* +1 dx <C{T h , f)S E l»*V»l(*T,«») / l|£>xV h ||dx 

JT n=0 VT€ Q ™ 

* /■ 

= C(T h ,f)^ / l'*V>l(*T.*n)ll^»V fc || rfx . 

n=0 VT€Q JT 


n- 0 VTGQ 


(4.20) 


The function |v?| is bounded and continuous on a bounded domain so in the limsup^o limit, the right- 
hand-side integral in Eqn. (4.20) approaches the measure-valued function (\<p\,p) which completes Lemma 
4.3. □ 

LEMMA 4.4. Let 6 Co(R d x R+) and assume that v* satisfies the conditions of Theorem 3.L 

The following measure-valued bound exists for h — > 0 


lim sup 
*-> o h 


51 A 4 E 7r *^T,<n) E / 5^(*,t)fW*,*)) dxdt 

n = 0 VT£Q Mi€T */R d xR+ i=1 1 


< 


C(T fc ,f) 
(fc + 1)! 




(4.21) 


where C{Th,$) is a constant that depends on Th and ||Z?* +1 f v ||. 

Proof. Choose c p(x,t) such that supp(v?) C Q x [0,r]. Recall that represents an 

approximation of the flux divergence integrated in a simplex T. By direct calculation 


N /d + 1 

TThvKzTi t n ) | 5^ $" T + e 

n=0 VT€Q 


U=1 


\ N , d 

-5>E / j ^n) ^ ^ ^ x» 

/ n= 0 VTSQ*' T »=1 

iV - <f 

= E At E / 7r '*^( X - t n)ECx i ( V '*) dx 

n = o vreQ*' 7 ' *-i 

N r d 

+ E A< E / _7r hV 7 ( x > f n))X] f ’,*i( V; ‘) 

n=0 VTeQ'' T i=l 


(4.22) 


dx 

(4.23) 


where is the quadrature error in calculating the flux divergence. From Assumption (H2), this quadrature 
error is assumed to be of the form 


lki* ) ll = 


r d d+1 

,/T i=l i=l 


< C(T h ) 


h k+i 

(* + !)! 


D k z +l 


K)K), Xj 


consequently for i r/»<p(x, £) bounded by a constant absorbed into C{Th) 


N 


E a * E TT h <p(XT,tn)t 

n=0 VTeQ 


N 


h k+1 ^ 

<C(Th) 


(t+1)! 


n=0 


r EAiWIW). 


u=l 


(4.24) 


(4.25) 
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Combining this result with Lemma 4.3 formally bounds the quadrature error term 


lim sup 

/i->0 


N 


n—0 VTeQ 




(4.26) 


(* + l)! 

Next, apply Green’s formula in each simplex to the first right-hand-side sum appearing in Eqn. (4.23) 

dx dt (4.27) 


E At E f v k<p(x,t n )'J2t^ l (vh) dx = f E / 

n= 0 VTeQ*' 7 ' i=l n=0 * //n VTeQ * /T <=1 

= -E/.E / T Etf f,(v ‘>^‘ 

n=0 * // VTeQ JI *= 1 

n=0 ^ /n VTeQ t=l 


Recall that 717 ^ and f are both bounded and continuous functions. Upon utilizing Remark 1 (i) and the 
compact support of it follows that the second right-hand-side sum of Eqn. (4.28) vanishes identically. 
Examining the remaining right-hand-side term in Eqn. (4.28), observe that 


E/.E jC(£^(n>-tg«*M 


dx dt\ 


^ C °Y [ E [ n f ( v /») - f ( v )ii <*** 

n=0‘ ,/n VT€Q*' T 

4 - CiV [ Y [ \\DxVh<P - Dz<p\\ ||f(vOII dxdt 

n=0 J I\TeQ JT 

(4.29) 


The first right-hand-side sum of Eqn. (4.29) is equal to H^v*) - f( v )llL l (Qx[o,r]) converges to 0 as 
h — ► 0, see Remark 1 (iii). Since stays bounded and f continuous, f(v^) stays bounded by a constant. 
The second right-hand-side sum of Eqn. (4.29) is bounded from above by \\D x 7Tk^p — Ae^llLMGxjo.r]) w ^ich 
also converges to 0 as h — > 0. Thus, it is concluded that 


lim 

h->0 


E/.E 

n=0 J 1 VTeQ JI \i-l t=l / 


dx dt 


= 0 


and consequently 


l™E At E / T ^(x,*n)E4(v/.)^ = -E/„E j T Y^^ dxdt 

n=Q VT£Q JT i= 1 n =0' // VTeQ*^ *=1 


(4.30) 


(4.31) 


Considering the second right-hand-side sum term in Eqn. (4.23), from Remark 1 (ii) it follows that 


N 

\T At T 


[ n= 0 


vreQ 


, d 
/ (^(^(xr^n) - 7T^^(x,t n ))E^«) dx 
Jt i= 1 


< E 


where 


s = c E il E j 

r» wt r\ J l 


I *h<P(XT, t n ) ~ 7Qi<l(s,t n ) 


n=0 VT€Q 

Since ||D x 7 r^y?|| is assumed bounded by a constant, 
1 7r/,<p(xr, in) - * h<p(x,tn) 


h 


E Ivt-vti*. 

A/i,Af>€T 


Sr 


dx 


Sr 


D x (n h ip) ■ 


XT “ Z 


dx < C h d 


(4.32) 
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where C is independant of h. Inserting this bound yields 
N . d 

E A * E / \^h<fi(x T ,t n ) - n h <p{x,t n )\ 5^r xi (v2) dx 

n = o vt€Q‘' t »=i 


N 


<aT‘'I E IIVj-v.I* 

n= 0 VTeQMi,Xfj€T 

(4.33) 


so that 


lim 

/i-M) 


N p d 

E A * E / ( 7r W( X T, t n) ~ *hV( x ,t n )) 

n=0 VT6Q‘ /T »=1 


dx 


= 0 


(4.34) 


Rearrangement of the bounded terms as h -4 0 completes Lemma 4.4. □ 

Proof. (Theorem 3.1) Multiply Eqn. (2.9) by <p(M u t n ) |Cj|, where <p(x,t) is a test function in 
Cq (R d x [0, +oo[), such that supp(<p) C Q x [0, r]. Summation on the indices n and i over time slabs 
and vertices, respectively, yields 

E E |C<MM i( t n ) (W^-WH+^At E v(Mi,tn)*?.T^0 . (4.35) 

n= 0 MitTh n=0 Mi€T h T,Mi€T 

From Lemma 4.2, 

limV y \Ci\ifi(M it t n ) (W** +1 - Wj 1 ) = - [ ^w(v fc (x,t)) dx- [ ^(x,0)w(v o (x)) dx . 

/l_>0 n^o h'xR+dt J R , 

The space term is rewritten 

N N 

E A * E E = E Af E E 

n=0 Afj T;Af , €T n-0 A/iGT 

JV 

+ E A *E E MAfi.tn) - 7r h v3(x T ,t n ))$" T (4.36) 

n= 0 VT€QMi6T 

where once again 7r^(/?(xT, in) denotes the midpoint value of the linearly interpolated function for constant 
t n . Examining the first right-hand-side sum of Eqn. 4.36, recall the result of Lemma 4.4, 


lim sup 
*-►0 h 


(E A * E E *M*T,t n )*Z T + [ E J^(x.*)f i (v*(x,t)) dxdt\ < 

ln=0 VTcQA/iCT JR J xR+ j=1 x * J 


C(%, f) 
(* + l)! 




Next, examine the second right-hand-side sum of Eqn. 4.36. From boundedness of ||I7y?j| combined with 
Assumption (HI) 

E ll v ”- v ”ll ( 4 - 37) 

Mj €T 


yielding 


N 


E At E E 

n=0 VTeQMiGT 


<ChE A « E E ll # r,ill 

n=0 VT6QMi€T 

< ch d jrAt y y iiv;-v?h . ( 4 . 38 ) 

n=0 VT€QMi,MjeT 


Consequently, from Lemma 4.1 as h — ► 0, 

N 

E At E E ’ *») - *hV(XT, *n)) *" T 

n =o vreQM.er 
which completes the proof of Theorem 3.1. 0 


(4.39) 
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5. The “N” Residual Distribution Scheme. An important example of a residual distribution 
scheme is the “N” (Narrow) scheme. It was first considered by Roe [27, 28] and Deconinck [8] for scalar 
equations. Here we consider the system extension due to van der Weide and Deconinck [31] and generalize 
their scheme to symmetrizable conservation laws 


d 

w(v), t + E f< ( W ( V )).*i = 0 ' 
1=1 


(5.1) 


Repeating Eqn. (2.18) of Sect. 2.2, our general strategy is to utilize a piecewise linear representation of the 
entropy variables themselves so that v/, € V/,. In a simplex T 


UP< w ( v h)),x;^ dxz= E<^<) J T ( v h), Xi dx = \T\^2(Ai)T(yh), Xi ,\ 


(5.2) 


with 


NQ 

{ Ai ) = E w ' Ai(v h (qi )) ,qt eT (5.3) 

/=1 

computed using iVQj>oint numerical quadrature. For purposes of analysis, it is convenient to define the 
symmetric matrices Kj t r € R mxm 


Kj,T = <r(^)T , VMj G T 


(5.4) 


where n J T € are the inward pointing normal vectors of the face of simplex T opposite vertex Mj scaled 
by the integral measure of the face, see for example Fig. 4.1. Also define K ± = (K ± \K\)/2 where \K | is 
calculated in the usual matrix sense via eigensystem decomposition. Due to the scaling of vector normals, 
UvA/jeT^r = Consequently, we have that ]C VJV/ ^ er Kjj — 0 and the identity 

E Kfr = - E K7.T ■ < 5 ' 5 > 

VAf^eT VA/jGT 


For the set of matrices {A*} equal to the Jacobian matrices of the Euler equations evaluated at a single 
state, it is shown in [1] that (X)vA/,er ) ‘ s nonsingular everywhere except when the state corresponds to 
a stagnation point. More generally, if we define (formally) the matrix N E R mxm in a simplex T 

n t = { E > (5 6) 

VvA^er / 

it is shown in the same paper that the matrix product KjN , VMj € T and its appearance in the N-scheme 
always has meaning, even at stagnation points. Hence forward, we assume that the matrix Nt always exists 
in the sense just described. 

Using these definitions, one can easily derive the following relationship for $x 
d \ d _ 

dx = |T| E^Mv/O.zi.ir “ E Kj'TVj • (5.7) 

:= 1 / «=1 VM,€T 

Fundamental to the N-scheme is following decomposition formula for 

*i,T = K+ T (Vj - Vi?"™) , VMj G T 



(5.8) 
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which is often called “upwind” because it represents a generalization to R* 4 of two-point upwind differencing 
for model scalar advection. Perhaps surprisingly, the requirement that &j,T represent a decomposition of 
*r, i.e. 


*r = E *j,t = E > ( 59 ) 

VM;6T VA^GT 

uniquely determines Vj? flow when Nt exists 

yrnflow = _ Nt K~ T Vi . (5.10) 

VAfi€T 

After some rearrangement, the N-scheme decomposition formula can be written in the following compact 
form 


*,,T = Kt Nt E KJt (Vj - v { ), VM, e T . 

va/^gt 


(5.11) 


The N-scheme then evolves the solution in time using the algorithm given earlier by Eqn. 2.9. Repeating 
this algorithm again while taking care to indicate the underlying dependence on v ^ and the nodal degrees 
of freedom V that describe v^: 


N-scheme in Symmetrization Variables : For all Mi E Thy n > 0, and v/ 4 E Vh 

r = w r - ^ r e t h*b . v»« = *wr») (5 . 12) 

1 W? = w(v„(M,)) 

The primary interest in the N-scheme for approximating conservation laws centers about a local discrete 
maximum principle exhibited by the N-scheme for scalar advection equations in R d . To see this, let v/,, Vi, 
and Wi denote the scalar (m = 1) forms of v^, V*, and W* respectively. Consider a numerical solution at 
steady state = v £. From Eqn. 5.11, the nodal degree of freedom at vertex Mi satisfies 


o 

II 

M 

E -k+ t n t k- t (v; -v;) 

(5.13) 

VTGT h ;A/iGT Mj GT \Mj^Mi 


= £ 

E {Vi - Vj) , a' 3 > 0 . 

(5.14) 


VTeT^Af.ET Mj 


This latter equation implies a local discrete maximum principle. More precisely, let adj a {M{) denote the set 
of vertices adjacent to Mi with nonzero weights a, then VMj E Th 


min V* < V* < max V.* . 

Mj Gadj a ( A/i) J Mj6*d} a (M t ) J 


Examining the time dependent problem in the scalar (m = 1) case, one easily derives a similar maximum 
principle result for n > 0 


min (V7\V7*) < vy* +1 < max (V. n ,V?) 

Mj €adj a ( M, ) 3 - A^eadj^M;) 3 


under the CFL-like condition at each t n 


At < max 


\C\i 


V Mi G Th X^VTgTmA/.GT Sa/jGTjA/j^^ -k+ t n t k- t 
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6. Energy and Entropy Analysis of the System N-Scheme. In this section, an energy analysis 
of the system (m > 1) N-scheme is undertaken as first described in Barth [5]. As a motivational exercise, we 
first consider analysis of the linear (constant coefficient) form of Eqn. (2.4). We then analyze the energy and 
entropy consequences of the N-scheme for nonlinear systems of conservation laws. This latter analysis shows 
that the N-scheme using symmetrization variables and numerical quadrature satisfies an entropy inequality 
in the limit of mesh refinement. 

6.1. Energy Analysis of the System N-Scheme: the Linear Case. In the linear system case, 
the numerical scheme (5.12) can be viewed abstractly as an Euler explicit integration of the semi-discrete 
matrix equation 

DV yt +LV = 0 (6.1) 

where V £ W is a vector representing the nodal degrees of freedom, D £ W xs a symmetric positive definite 
(SPD) matrix, and L £ R sxs a general real-valued matrix. The energy evolution is then given by 

i(V T Z?V),t + V r LV = 0 , L = (L + L T )/2 . (6.2) 

where L denotes the symmetric part of L. Energy boundedness is demonstrated if it can be shown that the 
symmetric part of L is positive semi-definite, i.e. for all V 

V T ZV = V r LV >0 . (6.3) 


Now suppose that this abstract matrix equation originates from a discretization procedure such as the 
N-scheme. The total energy associated with the matrix L can be computed and assembled on an element- 
by-element basis 

V T Zv= V£Z t V t ( 6 . 4 ) 

T€T h 

where Vj and Lt denote the nodal degrees of freedom and element matrix associated with a simplex 
T. Consequently, to demonstrate energy boundedness of the abstract linear system it is sufficient but not 
necessary to show 

v£Z r Vt > o , VT £ % . (6.5) 

Turning attention now to the N-scheme. For ease of exposition, we will show the development in two 
space dimensions but the generalization to R d will be clear. Next, consider the linear (constant coefficient) 
form of Eqn. (2.4). In this linear model, the conservation and symmetrization variables are related by the 
constant matrix Ao, i.e. 

W i = A 0 V j , VMj £ Th ■ (6-6) 


The SPD matrix D appearing in (6.1) wx>uld then be block diagonal with m x m blocks corresponding to 
each vertex Mi of the form |C|jAo. In two space dimensions, the system N-scheme decomposition (5.11) 
reduces to the following space discretization for a simplex T with local numbering T(M\, M 2 , A/3) 


LtV t — 



K+ 


Kt 


- 

K t 


*r 

T~ 

/ v , \ 

-F 


[TV] 

*2~ 


v 2 




K s’ 


V v 3 ) 


(6.7) 


with K ± symmetric and [N] a block diagonal matrix [N] — diag (N,N f N). The symmetric part of L is given 
by 


Kt 


Lt ~ 




*3 + 


1 

+ 2 


k r 


*r 

T 

1 

*r 

K+ 

[TV] 

k; 

+ 2 




£ 3 " 


£ 3 - 


[N] 


Kt 

Kt 

Kt 


( 6 . 8 ) 
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Examining rows of Lt or Lt , observe that the row sum is nonzero. However, we can add the following block 
diagonal matrix to the element matrix L 


1 

2 


K x 

K 2 

K 3 


(6.9) 


so that rows and columns of the Lt sum to zero. These additional terms have no impact on the constant 
coefficient discretization of the Cauchy problem. These terms all vanish identically when summed for all 
elements sharing a mesh vertex since the geometry surrounding the vertex is closed. Hence forward, we will 
include these terms in our definition of Lt and Lt yielding 


\K\i 

1 

+ 2 

Kf 


Kf 

T 

1 

+ 2 

Kf 


Kf 

\K\2 

Kf 

[N] 

Kf 

Kf 

[iV] 

Kf 

\k\ 3 

4 

Kf 


Kf 


Kf 


Kf 




Next, rewrite off-diagonal term such as 


Kf NKJ + KjNKf 


in the following form 


Kf NKJ + k~NK+ = kiNKj - kf N kf - kj N k~ 
Consequently, Lt can be rewritten as 


Lt= 2 


1 

+ 2 


k, 


K i 

T 


K-2 

[N] 

k 2 



. ** . 





' Kf 




Kf 

r 


kf 


- 

Kf 

[N] 



W + 


Kf 



Kf 

kf 

Kf 


1 

+ 2 

-Kf 

-Kf 



-Kf 

~Kf 

[N] 

~Kf 

-Kf 


-Kf 


~Kf 


~Kf 


( 6 . 10 ) 


( 6 . 11 ) 


Note that the first term appearing on the right hand side of Eqn. (6.11) gives rise to a quadratic form 
with positive energy so our only concern is the remaining terms on the right hand side on this equation. 
Before proving positive semi-definiteness of (6.11), we first review a simple result concerning the spectra of 
noncommuting matrices. 

LEMMA 6.1. The nonzero parts of the spectrum of AB and BA are identical for all matrices A G K mxn 
andBtW 1 *™. 

Proof. Omitted, see for example Axelsson [3] p. 69. □ 


Next we prove positive semi-definiteness of a specialized matrix in product form. 
Lemma 6.2 (Golub[12]). The matrix 


‘ A 

0 

0 ' 


‘ A ‘ 


' A ' 

0 

B 

0 

— 

B 

N 

B 

0 

0 

C 


C 


C 


N = [A + B + C\~ l 


is positive semidefinite for all A,B,C G K nxn symmetric positive definite. 
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Proof. Let 

and congruence transform L 

z _ 1 / 2 z,z _1/2 = 


z = 


.4 0 0 
0 £ 0 
0 0 C 


I n 


In 


In 


”1 

' A l/2 ' 


■ A 1 / 2 ' 

— 

£1/2 

N 

£!/2 


c 1 / 2 


. c 1 ' 2 . 


-i T 


= hn ~ P- 


Next use Lemma 6.1 concerning the spectra of nonsquare matrix products. In the present case Lemma 6.1 
this implies that 


Eigenvalues 


A 1 / 2 

B 1 / 2 

C x > 2 


N 


A l/2 

£1/2 

C 1 / 2 


= Eigenvalues^iV 1,/2 (.4 + £ + C)N 1//2 ^ + 2 n zeros 

= Eigenvalues (^N(A 4 B + C 4* 2 n zeros 
= Eigenvalues(/ n ) 4 2n zeros 


( 6 . 12 ) 


and consequently 


hn~P 


is positive semidefinite. From this result it follows immediately that 

L = Z 1 / 2 (/ 3n -P)Z 1 / 2 

is also positive semidefinite. □ 

The extension to A, B, C > 0 and (A + B + C) > 0 follows by considering the perturbed matrices A € = A+el, 
B t = B 4 e/, and C e = C 4 eJ and letting e | 0. 

Returning to the system N-scheme ? we now can prove the main result of this section. 

THEOREM 6.3. The system N-scheme discretization of the constant coefficient form of (2.4) is energy 
bounded with the element energy matrix (6.11) positive semi- definite, i.e. V T LV > 0. 

Proof. Since N = [K+ 4 K } 4 Kf]~ l = [—AT - A 2 “ - the result follows immediately after 

application of the Golub lemma to (6.11) together with Eqn. (6.4). □ 

6.2. Energy and Entropy Analysis of the System N-Scheme: the Nonlinear Case. In this 
section, an energy analysis of the N-scheme is presented for nonlinear systems of conservation laws. This 
energy also represents an approximation to the entropy inequality equation (2.2), see Hughes [15] or Barth 
[6, 4] for related entropy analysis of finite element discretizations. Specifically, we show convergence to 
an entropy inequality for the N-scheme with exact integration. We then show that with sufficient order 
numerical quadrature that the entropy inequality is retained in the limit of mesh refinement. 

Lemma 6.4. Under the assumptions of Theorem 3.1, the limit v of v* defined by the conservative system 
N-scheme satisfies 

d 

(6.13) 


± [ H(v)+ [ £<?>)* dS< 0 
at Jn Jd n 


Proof. Consider the system N-scheme decomposition (6.7). Unlike the constant coefficient linear case, 
the diagonal term (6.9) 


1 

2 


K i 


K 2 


k 3 
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can not be added to the element matrix Lt in the nonlinear case without changing the discretization. 
Consequently, the energy associated with the simplex T(M \, . . . , A/</+i) must include this term, i.e. 


1 d+i 

V?irVT = :J] V ^ V i + ( 3(V‘ V w) 


(6.14) 


i= 1 


where the quadratic form Q is positive by Theorem 6.3. The task is to show that in the limit h — ► 0 that 
the first right-hand-side term appearing in Eqn. (6.14) converges to 


r d 

/ VG'(v) Hi dS , 
Jdn i=1 


(6.15) 


the integral of the entropy flux. Recall that V describes the nodal degrees of freedom in the piecewise linear 
space Vh E V*. In a simplex T we have 


rf+ 1 rf+i _ 

Evj^v, = E ( v J^ v i - v ^ v 0 

i=i 


j=2 
d+1 

j=2 

d+1 

= E(V>+ v i) T *i (**»■&) 

J= 2 


(6.16) 

(6.17) 

(6.18) 


where /j* denotes the vector from vertex Af* to Mj. Thus we can define (by identification) the vectors 
Pj € R m such that 

d -\- 1 d r \ 

.§g. 


i=i 


i-i 


that is 


d+l 


p J = E z u( v ' +v ‘) r ^ 


1=2 


where we P X j is the j-th component of lij. Consequently, we obtain 

5 ivT*n = J T < (g^< v ‘)^) Lp- 


dx + er 


where G is the entropy flux associated with v and 

d 


er 


■/,( 


JjPl'a*) 


T 

~ Vfc 


j=l J 


da; . 


(6.19) 


(6.20) 


Hence, er can be estimated, 


j = 1 

d r 


^ - v h Aj( v ft) 


= < 




- v h ^it^) 


da: 


dx 


1/2 

• (/. 


1/2 

► /. 11^*11 


-> 1/2 

IPv,|| 2 dx 


J 


dx 
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because Dv/, is constant in each simplex T. Proceeding as in Lemma 4.1, we see that the function w defined 
by (see Figure 6.1 for a 2D illustration) 

t»|T= E (Vi + Vj) XD(J 

[i,j] edge of T 

where \d is the characteristic function of the set D converges in Lf oc to 2v when /i — > 0. By definition of 

D\i 2 



Fig. 6.1. Geometrical elements for the definition of w 


Kj , we have 


thus we see that 


cf-f" 1 d 

E^(v>-v 1 ) = m E^a^’ 

i=i 


i=2 


d /d + 1 


E(E'W)£Hn£^. 

i=i J 


(6.21) 

, *■' ’ / cfx; * — ' ' ax^ 

j=l \f=2 

Due to the boundedness of v^, we can apply the dominated convergence theorem and equation (6.21) thus 
yielding 


lMm 


1/2 


- v h 


dx 


for any bounded domain Q C K d . This ends the proof. □ 

Theorem 6 .5. Under the assumptions of Theorem 3.1 , the limit v o/v* defined by the system N-scheme 
satisfies the entropy enequality 

f ^H(x)dx+ [ E 7 r- Gi ( v ) dx + f <p(x,Q)v{x,0)dx < 0 . (6.22) 

Jn Ot J RJxR+^^i •/r<‘ 


/or <mj/ smooi/i > 0. 

Proof. The first observation is 


d ff(v) 
dt 


dw 

~dt 


(6.23) 


The second observation is that in a simplex T(Mi, M2, A/3) 

d+1 d+1 d+1 


Y.<PiV7*i ~ 


i=l 


1=1 


i=l 


(6.24) 
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where pc = Then consequently, 


d+i 1 d+l 

E V J^ = oE v M + ( 3 ( v *’- ’ Vd +^ - 


(6.25) 


j = i j=i 

where Q is positive. Thus we get, because Pj,pg > 0 

d+1 . d+l _ d+1 

^ViVjQi > tpc-YlvjKjVj + E^ -VoWlii ■ 


(6.26) 


1=1 


J = 1 


i=l 


The last observation is that in a simplex T 


d+l 

i=i 


d+i 


<C*IWI«, E H V *H Ill'll H v <- v jH 


(6.27) 


»,j= i 


where C is a constant depending on the mesh only and Cij = Kf N K } . From this we conclude that 


d+i 

Efo-Vo)Vf*i 


1=1 


< C'/i 2 HU^IIcx) 1 1 V, — Vj || 


(6.28) 


where C* depends on maxt||v/i|| which is uniformly bounded by assumption. Since the mesh is regular, 
h 2 < C n \T\ for a well chose constant independant of the mesh. Lemma 4.1 shows that B — » 0 when h — ► 0. 
Using the same arguments as in Theorem 3.1 and Lemma 6.4, we conclude that the semidiscrete scheme 
satisfies an entropy inequality. □ Combining the previous results of this section together, we finally conclude 
with the following Corollary: 

Corollary 6.6. Under the assumptions of theorem 3.1 and 6.5, the system N-scheme associated with 
the quadrature formula (2.20) satisfies in the limit h — > 0 an entropy inequality for any smooth ip > 0 


/ ^H(vMx + [ ^^ff(v)dx + / 

J Q ut JU d xU+ u x % Jr 


<e(x,0)v(x,0) dx < C-r— — y(M,/x) • ( 6 - 29 ) 

Jr<* + 1J- 


7. Numerical Results. In this section, numerical validation of Theorem 3.1 is provided via N-scheme 
calculation of smooth and discontinuous solutions of a scalar conservation law and system Euler equations 
for subsonic, transonic, and supersonic blunt body flows. Recall that Theorem 3.1 states, under classical 
assumptions, that numerical solutions of the N-scheme with adaptive quadrature converge to a function for 
which the residual 

f ( £) + y; 0 > ,fj(w(j, t)) ] dxdt+ f <j>(x,0)w o (x)dx 

Jqx[o,t] \ i=1 ) Jq 

may not vanish as in the classical Lax-Wendroff theorem. Instead, the residual is bounded by a measure- 
valued function with strength independent of the mesh size h such that the measure strength can be made 
arbitrarily small by making the number of quadrature points sufficiently large. As a practical matter, 
as will be shown in Sect. 7.1, the convergence is very rapid when derivatives of the flux components, 
f*, are well behaved. In addition, an adaptive quadrature scheme is proposed and tested which greatly 
reduces the computational cost of the N-scheme with quadrature. The adaptive quadrature strategy uses a 
simple estimate of solution smoothness to select the number of quadrature points thus producing an overall 
economical discretization method since most elements need only use one interior quadrature point (even for 
second order accurate extensions [1]). 
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7.1. ID Conservation Law. Consider the scalar Cauchy problem (1.1) 

f + (/(«),* = 0 for (x,t) G R x R+ 

\ u(x,0) = u 0 (x) 

First observe that the upwind scheme (1.5) of Sect. 1 on a uniform mesh can be rewritten as 

AXj + ($+—1/2 + *2+1/2) = 0 

with 

*7+1/2 - ( a )7+l/2( U J + 1 ~ U i ) > *7+1/2 = ^ a )7+l/2( U J +1 ~ U i) ( 7 ’ 2 ) 

and from Sect. 1 

(a), +i/2 s f(Uj+l) - /( - ^ = I' a(i ru(0)) ^ , ™(0 =«,-+* («j+i - «j) • (7-3) 

Note that on a nonuniform mesh, A xj is replaced by the lumped average A Xj = (Axj_i/ 2 + 1 / 2 )/% 

although other possible definitions are possible, e.g. Axj = (pt_^ 2 Axj_!/ 2 + Pj+i/ 2 ^ x j+ 1 / 2)1 Pj T 1/2 - 
(1 ± sgn((a)j^i/ 2 ))/2. Consistent with the previous analysis, our first numerical experiment implements a 
variant of this residual distribution upwind scheme of the form 

(*!-./>■ + *r + ./0=° <7 - 4) 

with residual distribution calculated via numerical quadrature 

NQ NQ 

$ 7+ i /2 = E w » a ( ,ru ^"^+ 1 *7 h /2 = a ( 7ru (q»)) + ( M j+i ~ u j ) ■ ( 7 - 5 ) 

/=i «=i 

The scheme (7.4) would be conservative if 

*7+1/2 + *7+1/2 = /(«;+») - /(«*) 

but due to the use of numerical quadrature 

NQ ,1 

*J +1/2 + *7+1/2 - ^2 Ui a (™(qi)) (u j+ i - Uj) ± / a(mi(0)d£{uj+i -Uj) = f(u j+ i) - f{uj) . 

1=1 Jo 

Even so, from Theorem 3.1 we still expect convergence to weak solutions provided sufficient order numerical 
quadrature is employed. 

7.1.1. ID Numerical Experiment: Fixed Gauss Quadrature on Nonuniform Mesh. We first 
test the scheme (7.4) with Euler explicit time advancement for the smooth flux formula and initial data 

f(u) = e u , ti 0 (x) = sin(27rx) 

on successively refined meshes (Ax = 10“ 2 , (1/2) 10“ 2 , (1/2) 2 10 -2 , (1/2) 3 10“ 2 and (1/2) 4 10 3 ). To 
eliminate superconvergent behavior of measured error norms due to mesh uniformity, we make the spacing 
between successive mesh points alternate between the values Ax and Ax/2. In evaluating the distribution 
formulas (7.5), iVQ-point Gauss quadrature formulas are used with 1 < NQ < 3 to validate Theorem 3.1. 
Selected results are given in Table 7.1 which tabulates L 1 , L 2 and L°° norms of the difference between the 
numerical solution u” given by standard conservative scheme (7.1) and the nonconservative provided 
by the scheme (7.4) on meshes with decreasing Ax at time nAt = 0.5 (after the shockwave has appeared). 
Figures 7.1(a-d) graph the solutions before and after the occurence of the shockwave for the conservative 
scheme and the nonconservative schemes using 1, 2, 3, 4 and 5 point Gauss quadrature on a mesh containing 
100 unknowns. All the solutions are virtually indistinguishable before the occurrence of the shockwave. It 
is only after the solution becomes discontinuous that the importance of sufficient order Gauss quadrature 
is visually seen and multiple quadrature points needed. The tabulated results reveal two effects addressed 
by the theory: (1) the L l error eventually stagnates when h — ► 0 for fixed order quadrature and (2) the L 1 
error decreases very rapidly with increasing NQ. In fact, upon closer inspection, this error decreases much 
more quickly than (jp + 1)! typical of Gauss quadrature, see Figure 7.2. 
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mesh size, Ax 

L* (u nc t^c) 

L 2 {u nc u c ) 

L°°(u nc Uq ) 

#quad pts, NQ 

0.100 10" 1 

0.63531 10“ 2 

0.25617 10“ 1 

0.15162 

i 

0.500 10“ 2 

0.67850 IO" 2 

0.38770 10" 1 

0.35783 

i 

0.250 lO" 2 

0.70532 IO" 2 

0.55376 IO" 1 

0.67416 

1 

| 0.125 10“ 2 

0.72127 IO" 2 

0.73250 IO" 1 

0.10491 10* 

1 

0.100 io- 1 

0.12402 IO" 4 

0.50233 10- 4 

0.30796 10-* 

2 

0.500 10“ 2 

0.14468 IO" 4 

0.83082 10- 4 

0.74799 IO" 3 

2 

0.250 IO" 2 

0.15648 IO -4 

0.12657 IO" 3 

0.14325 10“ 2 

2 

0.125 10" 2 

0.16296 IO' 4 

0.18732 IO" 3 

0.33085 IO -2 

2 

0.100 10" 1 

0.28748 10“' 

0.42536 10“' 

0.23256 10-° 

3 

0.500 IO" 2 

0.27937 IO" 7 

0.44600 IO" 7 

0.35562 IO" 6 

3 

0.250 IO" 2 

0.27315 IO" 7 

0.48398 10~ 7 

0.49170 10~ 6 

3 

0.125 10~ 2 

0.27017 IO" 7 

0.57222 IO -7 

0.93983 IO -6 

3 


Table 7.1 

Numerical results for ID Cauchy problem. Numerical error between the conservative calculation u c and the nonconserva- 
tive calculation u nc using NQ point Gauss quadrature. 




(a) Solution before the shockwave. (b) Solution after the shockwave. 




(c) Solution before the shockwave, zoom. (d) Solution after the shockwave, zoom. 


Fig. 7.1. Numerical N-scheme solutions for (1.1) and tx° = sin(27ri/i). Solutions before the formation of a shockwave ((a) 
and (c)) and solutions after the formation of a shockwave ((b) and (d)). 
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7.1.2, ID Numerical Experiment: Adaptive Gauss Quadrature on Nonuniform Mesh. The 
results of the previous ID numerical experiment of Sect. 7.1.1 show very negligible sensitivity to the num- 
ber of quadrature points whenever the solution is smooth. This observation suggests the following simple 
adaptive quadrature scheme which uses a nondimensional measure of solution gradient to estimate solution 
smoothness: 

» If Ui+1 < JAxJZ, then the solution is smooth. Compute ^++ 1/2 and ^j>i /2 with 

max(|uj|, |uj+i |) 

NQmin point quadrature. 

• Else, compute 'ip'j+x/o and V 7 + 1/2 NQ m ax point quadrature. 

Repeating the calculations of Sect. 7.1.1, Table 7.2 tabulates the corresponding numerical results using 
the adaptive parameters NQ m in — 1 and NQ ma x = 2,3 at the time T = 0.5 (after the formation of the 
shockwave). Note that in these calculations, nearly all cells required only NQ m i n — 1 quadrature point 


mesh size, Ax 

(u nc u c ) 

L‘ 2 {u nc - u c ) 

L°°(u nc u c ) 

A" Qmin 

NQ mai 

0 . 100 10- 1 f 


0.11351 10 _a 

0.69521 10 -3 

1 

2 

0.500 10~ 2 

0.27998 10" 4 

0.13526 10“ 3 

0.12200 10 -2 

1 

2 

0.250 10- 2 

0.21424 10- 4 

0.16100 10 -3 

0.18236 lO -2 

1 

2 

0.125 10- 2 

0.18526 10“ 4 

0.20768 lO" 3 

0.36673 10" 2 

1 

2 

0.100 10- 1 

IligsMiliaiiM 

BEHEEM 

0.38681 10~ 3 

1 

3 

0.500 10“ 2 

0.13642 10' 4 

0.52260 10“ 4 

0.47082 10' 3 

1 

3 

0.250 10“ 2 

0.57902 10“ 5 

0.34398 10“ 4 

0.38955 10~ 3 

1 

3 

0.125 lO" 2 

0.22238 lO -5 

0.20259 10- 4 

0.35797 10' 3 

1 

3 


~ — “ “ Table 7.2 

Error between the conservative scheme and the adaptive quadrature scheme at t = 0.5 tism$ 1, 2 , or 3 point Gauss 
quadrature. 


with only 3-5 cells requiring NQ max number of quadrature points. This results in a notable savings in 
computational resources. These numerical results indicate that the quality of the solutions is comparible 
with those of Table 7.1 with some reduced accuracy that would be improved by a more stringent criteria for 
quadrature adaptation. We have not run this case with a second order upwind scheme, but we believe that 
the same strategy could be used since quadrature with NQ m i n = 1 points is second order accurate. Intact, 
it can be shown formally that to recover second order accuracy, the “exact total residual 2 *b 1/2 
needs only be recovered up to second order accuracy to have a second order accuracte scheme, see [1]. Finally, 
we note that other tests have been carried out, for example with the flux f(u) = exp(u 2 ) with similar results. 
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7.2. 2D Conservation Laws. Next, we present 2D solutions of the Euler equations of gasdynamics 
assuming a perfect gas relationship. Throughout the remaining 2D numerical experiments, standard quadra- 
ture formulas for simplicies are utilized with weights and quadrature point locations as given in Table 7.3, 
see also [29]. 


NQ 

Accuracy 

Barycentric coordinates, qi 

weights, lji 

1 

om 

(1/3, 1/3, 1/3) 

1/2 

3 

om 

(1/2, 1/2,0) 

1/3 

6 


(0.816847572980459,0.091576213509771,0.091576213509771) 

(0.108103018168070,0.445948490915965,0.445948490915965) 

0.109951743655322 

0.223381589678011 

7 

om 

(1/3, 1/3, 1/3) “ t 

(0.797426985353087, 0.101286507323456, 0.101286507323456) 
(0.470142064105115,0.470142064105115,0.059715871789770) 

0.225 

0.125939180544827 

0.132394152788506 

16 

om 

(0.057104196114518,0.065466994555014,0.877428809330468) 
(0.2 768430 13638124, 0.050210123211370, 0.672946863150506) 
(0.583590432368917, 0.028912084224389, 0.387497483406694) 
(0.860240135656220, 0.009703785126946, 0.130056079216834) 
(0.057104196114518,0.311164552244357,0.631731251641125) 
(0.276843013638124, 0.238648659731443, 0.484508326630433) 
(0.583590432368917,0.137419104134574,0.278990463496509) 
(0.860240135656220, 0.046122079906452, 0.093637784437328) 
(0.057104196114518,0.631731251641125,0.311164552244357) 
(0.276843013638124, 0.484508326630433, 0.238648659731443) 
(0.583590432368917, 0.278990463496509, 0.137419104134574) 
(0.860240135656220, 0.093637784437328, 0.046122079906452) 
(0.0571041961 14518, 0.877428809330468, 0.065466994555014) 
(0.276843013638124, 0.672946863150506, 0.050210123211370) 
(0.583590432368917, 0.387497483406694, 0.028912084224389) 
(0.860240135656220, 0.130056079216834, 0.009703785126946) 

0.047136736386776 

0.070776135796160 

0.045168098564740 

0.010846451821051 

0.088370177044746 

0.132688432214078 

0.084679449043492 

0.020334519128958 

0.088370177044746 

0.132688432214078 

0.084679449043492 

0.020334519128958 

0.047136736386776 

0.070776135796160 

0.045168098564740 

0.010846451821051 


Table 7.3 

Quadrature points and weights. The missing quadrature points are obtained by cyclic permutation as needed. 


The significant computational savings obtained by adaptive numerical quadrature in ID suggest using 
a similar strategy in higher space dimensions were the savings is even more dramatic. For any simplex T\ a 
criterion must be developed which determines if the numerical solution is locally smooth or not. For efficiency 
reasons, this decision should ideally be made from the information available in T only. Let Sj denotes the 
(physical) entropy at node M } and hr the maximum length of the edges of T . We have implemented the 
following heuristic criterion for use in the adaptive quadrature strategy: 


• If min Mi Mj^T 


Si 


1 > ^JlvrJh, then the solution is smooth. Compute the N-scheme decomposi- 


tion using NQmin point quadrature. 

• Else, compute the N-scheme decomposition using NQ max point quadrature. 


7.2.1. 2D Numerical Experiments: Euler Equations on Mesh Triangulations. We first study 
the effect of the loss of conservation and the influence of the number of quadrature points for the N-scheme 
with quadrature. To achieve this goal, we have selected three test cases that are simple yet representative of 
different flow regimes: a subsonic flow test case, a transonic flow case with mild shockwaves, and a supersonic 
flow case over a blunt body which produces a strong bow shockwave. The solutions are compared to those 
obtained by the reference conservative N-scheme using Z variables. Our intent is not to assess the accuracy 
of the solution with respect to a mesh-converged solution, but rather to see how the loss of conservation 
affects the structure of the solution compared with the reference solution on the same mesh. In particular, 
we qualitatively and quantitatively compare the overall structure of conservative and nonconservative N- 
scheme solutions by examining representative cross-sectional and/or boundary data plots. Additionally, 
we examine the behavior of numerical solutions with adaptive mesh refinement for the cases containing 
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discontinuities where exact discrete conservation is normally very important as it insures proper solution 
jump approximation. Ideally, it would be illuminating to perform uniform mesh refinement in evaluating 
the adaptive quadrature N-scheme as h -* 0. This was done in the ID calculations. Unfortunately, this is 
prohibitively expensive in 2D so we rely on multiple levels of adaptive mesh refinement to approximate the 
/i — ► 0 process. 

7.2.2. Subsonic Flow Case. This case is taken from Dervieux [9]. It is a flow over a cylinder with 
a Mach number at infinity of M oo — 0.38 computed on a relatively coarse mesh containing 3010 simplicial 
elements. Under these flow conditions, the flow remains subsonic and devoid of solution discontinuities. 
Figures 7.3(a-d) show Mach number isolines for N-scheme calculations using 1, 3, and 7 point quadrature as 
well as the N-scheme using the conservative Z variables. Figure 7.4(a) shows Mach number isolines for the 
N-scheme using the adaptive quadrature procedure described earlier with parameter values iVQ mm = 1 and 
NQmax = 3. Figure 7.4(b) provides a quantitative comparison of pressures on the surface of the cylinder 
using all the fixed and adaptive quadrature formulas as well as the conservative Z variable N-scheme. As 




(c) 7 point quadrature 




(d) conservative Z variables 


Fig. 7.3. (a-d) Mach number isolines for N-scheme calculations using fixed 1, 3 , and 7 point quadrature and the conser- 
vative Z variables for the subsonic cylinder problem , Moo = 0.38, on a simplicial mesh containing 3010 elements. 


expected, all calculations show no discernible differences. This results confirm our analysis if we assume 
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(a) 1-3 point adaptive quadrature 


(b) Pressure along top-bottom line of symmetry 


Fig. 7.4. (a) Mach number isolines for N -scheme calculations using 1-3 point adaptive quadrature and (b) the resulting 
pressure along top-bottom line of symmetry for the subsonic cylinder problem, Moo = 0.38, on a simplicial mesh containing 
3010 elements . 


that the support of the measure jj, is concentrated near discontinuities in the solution. Since there are no 
discontinuities in this flow and our quadrature formulas are at least second order accurate using single point 
quadrature, a Lax Wendroff theorem is satisfied up to 0(/i 2 ). 

7.2.3. Transonic Flow Case. The second 2D test case consists of transonic flow, Moo = 0.85, over 
the NACA0012 geometry with flow incidence of 1° computed on a baseline simplicial mesh containing 5050 
elements. The flow solution consists of both upper and lower surface shockwaves. Due to the 1° flow 




(a) Fixed 1,3, and 7 point quadrature and conservative (b) Three levels adaptive mesh refinement 
Z variables with 1-3 point adaptive quadrature 

Fig. 7.5. Surface pressure coefficient distribution on the NACA0012 airfoil using the N-scheme with (a) 1, 3, and 7 point 
quadrature and the conservative Z variable and (b) three levels of adaptive mesh refinement together with 1-3 point adaptive 
quadrature. 

incidence, the upper surface shockwave is notably stronger than the lower surface shockwave. N-scheme 
calculations were performs using fixed 1, 3, and 7 quadrature point formulas as well as the conservative 
Z variables on the baseline simplicial mesh containing 5050 elements. In addition, three levels of adaptive 
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mesh refinement were performed and solutions computed using the adaptive quadrature N-scheme with 
NQ m i n = 1 and NQ max — 3. Surface pressure coefficient values are graphed in Fig. 7.5(a) and Mach 
number isocontours shown in Figs. 7.6(a-d) for N-scheme calculations using 1, 3, and 7 point quadrature 
and the conservative Z variables on fixed mesh composed of 5050 elements. Similarly, surface pressure 
coefficient values are graphed in Fig. 7.5(b) and Mach number isocontours in Figs. 7.7(a-d) using three 
levels of shockwave-adapted mesh refinement together with the 1-3 adaptive quadrature point form of the 
N-scheme. Although the Mach number isocontour plots look visually very similar, the Fig. 7.5(a) graph 






FlG. 7.6. (a-d) Mach number isolines for N-scheme calculation using fixed 1, 3, and 7 point quadrature and the conservative 
Z variables for the transonic NACA001 2 problem , Moo = 0.85 and 1° flow incidence, on a simplicial mesh containing 5050 
elements. 

of the pressure coefficient on the body of the NACA0012 airfoil is more revealing. This graph shows that 
the location of the shockwaves depends on the number of quadrature points. Specifically, the use of single 
point quadrature leads to a significant change in shockwave location when compared to 3 and 7 point 
quadrature as well as the conservative Z variable scheme. For this particular flow, the effect of the measure 
H is not sufficiently reduced using one quadrature point but using three or more quadrature points seems 
sufficient to reduce conservation error less than truncation errors present in the conservative N-scheme. The 
comparability of the 3 and 7 point quadrature with the conservative Z variable scheme once again suggests 
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(a) 1 level adaptive mesh refinement, 
1-3 point adaptive quadrature 


(b) 2 levels adaptive mesh refinement, 
1-3 point adaptive quadrature 




(a) 3 levels adaptive mesh refinement, 
1-3 point adaptive quadrature 


(b) 3 levels adaptive mesh refinement, 
conservative Z variables 



Fig. 7.7. (a-d) Mach number isolines for N-scheme calculations using 1, 2, and 3 levels of adaptive mesh refinement 
and 1-3 point adaptive quadrature and the reference conservative Z variables on the 3 level refined mesh, for the transonic 
NACA0012 problemj Moo = 0.85 and 1° flow incidence. 


an adaptive quadrature implementation. The calculations presented in Fig. 7.5(b) are intended at checking 
whether the errors generated by the loss of conservation on refined meshes dominate the truncation error, 
even in an adaptive quadrature setting. With adaptive mesh refinement, all the computations in Fig. 7.5(b) 
are very comparable which further validates Theorem 3.1 and our adaptive quadrature strategy. 

7.2.4. Supersonic Blunt Body Flow. The last 2D test case consists of supersonic flow, Af©© = 3.5, 
over a circular cylinder geometry computed on a baseline simplicial mesh containing 4075 elements. The flow 
solution consists of strong bow shock forward of the cylinder geometry. Figures 7 .8 (a-d) show Mach number 
isocontours for numerical solutions computed using 1,3, and 7 point quadrature and conservative Z variable 
forms of the N-scheme. In addition, Mach number isocontours for 1-3 and 1-7 point adaptive quadrature 
N-scheme calculations are shown in Figs. 7.9 (a-b). Both fixed and adaptive quadrature calculations are 
compared in Fig. 7.9 for pressure data along the top-bottom line of symmetry. This latter figure shows a 
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(a) 1 point quadrature (b) 3 point quadrature 


(c) 7 point quadrature 


(d) conservative Z variables 





Fig. 7.8. Mach number isolines for N -scheme calculations using fixed 1, 3, and 7 point quadrature and the conservative 
Z variables for the supersonic cylinder problem, Moo = 3.5 on the baseline simplicial mesh containing 4075 elements. 



(a) 1-3 point adaptive 
quadrature 



(b) 1-7 point adaptive 
quadrature 



(c) Pressure along line of top-bottom symmetry 


Fig. 7.9. (u-b) Mach number isolines for N-scheme calculations using 1-3 and 1-7 point adaptive quadrature and (c) 
comparison of all fixed and adaptive quadrature N-scheme calculations along the top-bottom line of symmetry for the supersonic 
cylinder problem , Moo — 3.5 on the baseline simplicial mesh containing 4075 elements. 
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large difference between the one quadrature calculation and the other calculations. This difference is also 
clearly seen in the Mach number isocontour plot, Fig. 7.8(a). Perhaps more importantly, Fig. 7.9 shows that 
the 3 point quadrature (and 1-3 point adaptive quadrature) also produced incorrect shockwave locations on 
the baseline mesh, although the error is much smaller than that obtained using 1 point quadrature. Recall 
that for the transonic flow problem, 3 point quadrature was of sufficient order on the baseline and adaptively 
refined meshes for computing correct shock locations. Using 7 point fixed quadrature and 1-7 point adaptive 
quadrature yields solution shockwave positions in agreement with the conservative scheme. These results 
are also in agreement with the inequality (3.8) of Theorem 3.1, since the strength of the measure depends 
on the number of quadrature points but also on the supremum of a norm of higher derivatives of the flux, 
£)J +1 f v . Estimation of this norm is difficult, but it is reasonable that this number tends to infinity as the 
maximum Mach number also tends to oo. However, since the flux f is analytical in v and the Mach number 
finite, the right-hand-side of (3.8) still converges to zero, albeit more slowly. 

Next, we examine the effect of adaptive mesh refinement. Figures 7.10(a-c) show Mach number iso- 
contours for the N-scheme calculations using 16 point quadrature, 1-16 point adaptive quadrature, and 
conservative Z variables. Figure 7.11 shows a graph of Mach number along the top-bottom symmetry line 
for these same schemes as well as 1-7 point adaptive quadrature. Surprisingly, Fig. 7.11 shows small differ- 



FiG. 7.10. (o-b) Mach number isolines for N-scheme calculations using 16 point fixed and 1-16 point adaptive quadrature 
and (c) conservative Z variables for the supersonic cylinder problem , Moo = 3.5, on the baseline simplicial mesh with 3 levels 
of adaptive mesh refinement . 


ences in shock profile using 1-7 point adaptive quadrature for this problem with three levels of adaptive mesh 
refinement. It is only with 16 point fixed or adaptive quadrature that the adaptive N-scheme solutions match 
the conservative Z variable N-scheme. This demonstrates some slight dependency on the mesh parameter h 
not captured by the present analysis. 

8. Concluding Remarks. A number of upwind schemes are derived in quasilinear form and discrete 
conservation obtained by devising specialized mean- value linearized coefficients. This approach is problematic 
for systems such as magnetohydrodynamics, Euler equations with certain forms of chemistry, etc. where these 
specialized mean- values linearizations may not exit in closed form. In the present analysis, we consider a 
more general construction of these upwind schemes which avoids explicitly constructing these exact mean- 
value linearizations. Our construction is well-tailored to systems of conservation laws with convex entropy 
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Fig. 7.11. Comparison of the Mach number along a top-bottom symmetryline using the N-scheme with 1 
adaptive quadrature , 16 point quadrature , and the conservative Z variables. 


16 point 


extension. Using the tools of weak-* convergence, a Lax-Wendroff theorem has been derived for this class 
of nonconservative schemes utilizing numerical quadrature. By using sufficient order numerical quadrature, 
we show that correct weak solutions of conservation laws are obtained. Numerical results confirm the basic 
analysis but do show some weak interdependence of the mesh space parameter h and the required order of 
accuracy of the numerical quadrature. This indicates that further investigation and quantification of this 
effect is needed. 
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